      subroutine adaptw(nxp, nyp, nzp, i1, i2, j1, j2, k1, k2, toroid, strait, &
         rmaj, x, y, z, wgrid, w_soln) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp 
      integer , intent(in) :: nyp 
      integer  :: nzp 
      integer , intent(in) :: i1 
      integer , intent(in) :: i2 
      integer , intent(in) :: j1 
      integer , intent(in) :: j2 
      integer , intent(in) :: k1 
      integer , intent(in) :: k2 
      real(double) , intent(in) :: toroid 
      real(double) , intent(in) :: strait 
      real(double) , intent(in) :: rmaj 
      real(double) , intent(in) :: x(nxp,nyp,*) 
      real(double) , intent(in) :: y(nxp,nyp,*) 
      real(double) , intent(in) :: z(nxp,nyp,*) 
      real(double) , intent(inout) :: wgrid(nxp,nyp,*) 
      real(double) , intent(in) :: w_soln(nxp,nyp,*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, j, i 
      real(double) :: ws, rminor, rmajor, sinphi, cosphi, stheta, ctheta, phi, &
         theta 
!-----------------------------------------------
!
      do k = k1, k2 
         do j = j1, j2 
            do i = i1, i2 
!
               ws = toroid*sqrt(x(i,j,k)**2+y(i,j,k)**2) - rmaj + strait*x(i,j,&
                  k) 
               rminor = sqrt(ws**2 + z(i,j,k)**2) 
               rmajor = toroid*sqrt(x(i,j,k)**2+y(i,j,k)**2) + strait 
!
               sinphi = y(i,j,k)/rmajor 
               cosphi = x(i,j,k)/rmajor 
!
               stheta = ws/(rminor + 1.E-30) 
               ctheta = z(i,j,k)/(rminor + 1.E-30) 
               ctheta = amin1(ctheta,1.) 
               ctheta = amax1(ctheta,-1.) 
               cosphi = amin1(cosphi,1.) 
               cosphi = amax1(cosphi,-1.) 
!
!
               phi = acos(cosphi) 
               theta = acos(ctheta) 
!
!      wgrid(i,j,k)=wgrid(i,j,k)*(1.0+1.e3*cosphi**4)
!      wgrid(i,j,k)=wgrid(i,j,k)
!     &     *(1.+025.*exp(-100.*(rminor-0.5)**2))
!      wgrid(i,j,k)=wgrid(i,j,k)
!      wgrid(i,j,k)=wgrid(i,j,k)*(1.+0.75*cos(2.*theta-3.*phi))
!
               wgrid(i,j,k) = wgrid(i,j,k)*w_soln(i,j,k) 
!
            end do 
         end do 
      end do 
!
      return  
      end subroutine adaptw 


      subroutine contras(iwid, jwid, kwid, ncells, ijkcell, x, y, z, tsix, tsiy&
         , tsiz, etax, etay, etaz, nux, nuy, nuz) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(in) :: ncells 
      integer , intent(in) :: ijkcell(*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double) , intent(out) :: tsix(*) 
      real(double) , intent(out) :: tsiy(*) 
      real(double) , intent(out) :: tsiz(*) 
      real(double) , intent(out) :: etax(*) 
      real(double) , intent(out) :: etay(*) 
      real(double) , intent(out) :: etaz(*) 
      real(double) , intent(out) :: nux(*) 
      real(double) , intent(out) :: nuy(*) 
      real(double) , intent(out) :: nuz(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, ipjk, ijpk, imjk, ijmk, ijkp, ijkm 
      real(double) :: jacob, x1, x2, x3, y1, y2, y3, z1, z2, z3, rjacob 
!-----------------------------------------------
!
!     a routine to calculate the contravariant mesh
!     vectors at the vertices of the mesh
!
!
!
!
!dir$ ivdep
      do n = 1, ncells 
!
         ijk = ijkcell(n) 
         ipjk = ijk + iwid 
         ijpk = ijk + jwid 
         imjk = ijk - iwid 
         ijmk = ijk - jwid 
         ijkp = ijk + kwid 
         ijkm = ijk - kwid 
!
!     compute covariant base vectors
!
         x1 = 0.5*(x(ipjk)-x(imjk)) 
         x2 = 0.5*(x(ijpk)-x(ijmk)) 
         x3 = 0.5*(x(ijkp)-x(ijkm)) 
!
         y1 = 0.5*(y(ipjk)-y(imjk)) 
         y2 = 0.5*(y(ijpk)-y(ijmk)) 
         y3 = 0.5*(y(ijkp)-y(ijkm)) 
!
         z1 = 0.5*(z(ipjk)-z(imjk)) 
         z2 = 0.5*(z(ijpk)-z(ijmk)) 
         z3 = 0.5*(z(ijkp)-z(ijkm)) 
!
!     compute determinant of the metric tensor
!
         jacob = x1*(y2*z3 - y3*z2) + y1*(z2*x3 - z3*x2) + z1*(x2*y3 - x3*y2) 
!
         rjacob = 1./jacob 
!
!     compute contravariant base vectors
!
         tsix(ijk) = (y2*z3 - y3*z2)*rjacob 
         tsiy(ijk) = (z2*x3 - z3*x2)*rjacob 
         tsiz(ijk) = (x2*y3 - x3*y2)*rjacob 
!
         etax(ijk) = (y3*z1 - y1*z3)*rjacob 
         etay(ijk) = (z3*x1 - z1*x3)*rjacob 
         etaz(ijk) = (x3*y1 - x1*y3)*rjacob 
!
         nux(ijk) = (y1*z2 - y2*z1)*rjacob 
         nuy(ijk) = (z1*x2 - z2*x1)*rjacob 
         nuz(ijk) = (x1*y2 - x2*y1)*rjacob 
!
!
      end do 
!
!
      return  
      end subroutine contras 


      subroutine gridsolv(x, y, z, w, tsix, tsiy, tsiz, etax, etay, etaz, nux, &
         nuy, nuz, iwid, jwid, kwid, ijkcell, ncells, dx1, dx2, dx3) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(in) :: ncells 
      integer , intent(in) :: ijkcell(*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double) , intent(in) :: w(*) 
      real(double) , intent(in) :: tsix(*) 
      real(double) , intent(in) :: tsiy(*) 
      real(double) , intent(in) :: tsiz(*) 
      real(double) , intent(in) :: etax(*) 
      real(double) , intent(in) :: etay(*) 
      real(double) , intent(in) :: etaz(*) 
      real(double) , intent(in) :: nux(*) 
      real(double) , intent(in) :: nuy(*) 
      real(double) , intent(in) :: nuz(*) 
      real(double) , intent(out) :: dx1(*) 
      real(double) , intent(out) :: dx2(*) 
      real(double) , intent(out) :: dx3(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, ipjk, ipjpk, ijpk, imjpk, imjk, imjmk, ijmk, ipjmk, &
         ijkp, ipjkp, ipjpkp, ijpkp, imjpkp, imjkp, imjmkp, ijmkp, ipjmkp, ijkm&
         , ipjkm, ipjpkm, ijpkm, imjpkm, imjkm, imjmkm, ijmkm, ipjmkm 
      real(double) :: lama, g11, g12, g13, g22, g23, g33, wip, wim, wjp, wjm, &
         wkp, wkm, wavg, x11, x12, x13, x22, x23, x33, y11, y12, y13, y22, y23&
         , y33, z11, z12, z13, z22, z23, z33, rax, ray, raz, rg 
!-----------------------------------------------
!
!     a mesh generator based on
!     Winslow's variable diffusion
!     functional
!
!
!
!
!dir$ ivdep
      do n = 1, ncells 
!
         ijk = ijkcell(n) 
         ipjk = ijk + iwid 
         ipjpk = ijk + iwid + jwid 
         ijpk = ijk + jwid 
         imjpk = ijk - iwid + jwid 
         imjk = ijk - iwid 
         imjmk = ijk - iwid - jwid 
         ijmk = ijk - jwid 
         ipjmk = ijk + iwid - jwid 
!
         ijkp = ijk + kwid 
         ipjkp = ijk + iwid + kwid 
         ipjpkp = ijk + iwid + jwid + kwid 
         ijpkp = ijk + jwid + kwid 
         imjpkp = ijk - iwid + jwid + kwid 
         imjkp = ijk - iwid + kwid 
         imjmkp = ijk - iwid - jwid + kwid 
         ijmkp = ijk - jwid + kwid 
         ipjmkp = ijk + iwid - jwid + kwid 
!
         ijkm = ijk - kwid 
         ipjkm = ijk + iwid - kwid 
         ipjpkm = ijk + iwid + jwid - kwid 
         ijpkm = ijk + jwid - kwid 
         imjpkm = ijk - iwid + jwid - kwid 
         imjkm = ijk - iwid - kwid 
         imjmkm = ijk - iwid - jwid - kwid 
         ijmkm = ijk - jwid - kwid 
         ipjmkm = ijk + iwid - jwid - kwid 
!
!     compute covariant base vectors
!
!
!     calculate elements of the metric tensor
!
         g11 = tsix(ijk)**2 + tsiy(ijk)**2 + tsiz(ijk)**2 
         g12 = tsix(ijk)*etax(ijk) + tsiy(ijk)*etay(ijk) + tsiz(ijk)*etaz(ijk) 
         g13 = tsix(ijk)*nux(ijk) + tsiy(ijk)*nuy(ijk) + tsiz(ijk)*nuz(ijk) 
         g22 = etax(ijk)**2 + etay(ijk)**2 + etaz(ijk)**2 
         g23 = etax(ijk)*nux(ijk) + etay(ijk)*nuy(ijk) + etaz(ijk)*nuz(ijk) 
         g33 = nux(ijk)**2 + nuy(ijk)**2 + nuz(ijk)**2 
!
!     compute derivatives
!
         wip = 0.5*(w(ipjk)+w(ijk)) 
         wim = 0.5*(w(ijk)+w(imjk)) 
!
         wjp = 0.5*(w(ijpk)+w(ijk)) 
         wjm = 0.5*(w(ijk)+w(ijmk)) 
!
         wkp = 0.5*(w(ijkp)+w(ijk)) 
         wkm = 0.5*(w(ijk)+w(ijkm)) 
!
         wavg = (wip + wim + wjp + wjm + wkp + wkm)/6. 
!
         x11 = wip*(x(ipjk)-x(ijk)) - wim*(x(ijk)-x(imjk)) 
!
         x12 = 0.25*(wip*(x(ipjpk)-x(ipjmk)+x(ijpk)-x(ijmk))-wim*(x(ijpk)-x(&
            ijmk)+x(imjpk)-x(imjmk))+wjp*(x(ipjpk)-x(imjpk)+x(ipjk)-x(imjk))-&
            wjm*(x(ipjk)-x(imjk)+x(ipjmk)-x(imjmk))) 
!
         x13 = 0.25*(wip*(x(ipjkp)-x(ipjkm)+x(ijkp)-x(ijkm))-wim*(x(ijkp)-x(&
            ijkm)+x(imjkp)-x(imjkm))+wkp*(x(ipjkp)-x(imjkp)+x(ipjk)-x(imjk))-&
            wkm*(x(ipjk)-x(imjk)+x(ipjkm)-x(imjkm))) 
!
         x22 = wjp*(x(ijpk)-x(ijk)) - wjm*(x(ijk)-x(ijmk)) 
!
         x23 = 0.25*(wjp*(x(ijpkp)-x(ijpkm)+x(ijkp)-x(ijkm))-wjm*(x(ijkp)-x(&
            ijkm)+x(ijmkp)-x(ijmkm))+wkp*(x(ijpkp)-x(ijmkp)+x(ijpk)-x(ijmk))-&
            wkm*(x(ijpk)-x(ijmk)+x(ijpkm)-x(ijmkm))) 
!
         x33 = wkp*(x(ijkp)-x(ijk)) - wkm*(x(ijk)-x(ijkm)) 
!
         y11 = wip*(y(ipjk)-y(ijk)) - wim*(y(ijk)-y(imjk)) 
!
         y12 = 0.25*(wip*(y(ipjpk)-y(ipjmk)+y(ijpk)-y(ijmk))-wim*(y(ijpk)-y(&
            ijmk)+y(imjpk)-y(imjmk))+wjp*(y(ipjpk)-y(imjpk)+y(ipjk)-y(imjk))-&
            wjm*(y(ipjk)-y(imjk)+y(ipjmk)-y(imjmk))) 
!
         y13 = 0.25*(wip*(y(ipjkp)-y(ipjkm)+y(ijkp)-y(ijkm))-wim*(y(ijkp)-y(&
            ijkm)+y(imjkp)-y(imjkm))+wkp*(y(ipjkp)-y(imjkp)+y(ipjk)-y(imjk))-&
            wkm*(y(ipjk)-y(imjk)+y(ipjkm)-y(imjkm))) 
!
         y22 = wjp*(y(ijpk)-y(ijk)) - wjm*(y(ijk)-y(ijmk)) 
!
         y23 = 0.25*(wjp*(y(ijpkp)-y(ijpkm)+y(ijkp)-y(ijkm))-wjm*(y(ijkp)-y(&
            ijkm)+y(ijmkp)-y(ijmkm))+wkp*(y(ijpkp)-y(ijmkp)+y(ijpk)-y(ijmk))-&
            wkm*(y(ijpk)-y(ijmk)+y(ijpkm)-y(ijmkm))) 
!
         y33 = wkp*(y(ijkp)-y(ijk)) - wkm*(y(ijk)-y(ijkm)) 
!
         z11 = wip*(z(ipjk)-z(ijk)) - wim*(z(ijk)-z(imjk)) 
!
         z12 = 0.25*(wip*(z(ipjpk)-z(ipjmk)+z(ijpk)-z(ijmk))-wim*(z(ijpk)-z(&
            ijmk)+z(imjpk)-z(imjmk))+wjp*(z(ipjpk)-z(imjpk)+z(ipjk)-z(imjk))-&
            wjm*(z(ipjk)-z(imjk)+z(ipjmk)-z(imjmk))) 
!
         z13 = 0.25*(wip*(z(ipjkp)-z(ipjkm)+z(ijkp)-z(ijkm))-wim*(z(ijkp)-z(&
            ijkm)+z(imjkp)-z(imjkm))+wkp*(z(ipjkp)-z(imjkp)+z(ipjk)-z(imjk))-&
            wkm*(z(ipjk)-z(imjk)+z(ipjkm)-z(imjkm))) 
!
         z22 = wjp*(z(ijpk)-z(ijk)) - wjm*(z(ijk)-z(ijmk)) 
!
         z23 = 0.25*(wjp*(z(ijpkp)-z(ijpkm)+z(ijkp)-z(ijkm))-wjm*(z(ijkp)-z(&
            ijkm)+z(ijmkp)-z(ijmkm))+wkp*(z(ijpkp)-z(ijmkp)+z(ijpk)-z(ijmk))-&
            wkm*(z(ijpk)-z(ijmk)+z(ijpkm)-z(ijmkm))) 
!
         z33 = wkp*(z(ijkp)-z(ijk)) - wkm*(z(ijk)-z(ijkm)) 
!
!
         rax = g11*x11 + g12*x12 + g13*x13 + g22*x22 + g23*x23 + g33*x33 
!
         ray = g11*y11 + g12*y12 + g13*y13 + g22*y22 + g23*y23 + g33*y33 
!
         raz = g11*z11 + g12*z12 + g13*z13 + g22*z22 + g23*z23 + g33*z33 
!
         rg = -1./(2.*wavg*(g11 + g22 + g33)) 
         dx1(ijk) = -rax*rg 
         dx2(ijk) = -ray*rg 
         dx3(ijk) = -raz*rg 
!
!
      end do 
!
!
      return  
      end subroutine gridsolv 


      subroutine meshgen(nxp, nyp, nzp, x, y, z, w, w_soln, vol, dt, taug, &
         numit, ax, ay, az, bx, by, bz, cx, cy, cz, lama, lamb, lamc, dx1, dx2&
         , dx3, delta, psi, tsix, tsiy, tsiz, etax, etay, etaz, nux, nuy, nuz, &
         divaa1, divbb2, divcc3, iwid, jwid, kwid, ijkcell, istep, iota, &
         periodic, toroid, strait, rmaj, rwall, bzi, q0, cdlt, sdlt, dz) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nxp 
      integer  :: nyp 
      integer  :: nzp 
      integer , intent(in) :: numit 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      real(double)  :: dt 
      real(double)  :: taug 
      real(double)  :: lama 
      real(double)  :: lamb 
      real(double)  :: lamc 
      real(double)  :: delta 
      real(double)  :: psi 
      real(double)  :: toroid 
      real(double)  :: strait 
      real(double)  :: rmaj 
      real(double)  :: rwall 
      real(double)  :: bzi 
      real(double)  :: q0 
      real(double)  :: cdlt 
      real(double)  :: sdlt 
      real(double)  :: dz 
      logical , intent(in) :: periodic 
      integer  :: ijkcell(*) 
      integer  :: istep(*) 
      integer  :: iota(*) 
      real(double)  :: x(*) 
      real(double)  :: y(*) 
      real(double)  :: z(*) 
      real(double)  :: w(*) 
      real(double)  :: w_soln(*) 
      real(double)  :: vol(*) 
      real(double)  :: ax(*) 
      real(double)  :: ay(*) 
      real(double)  :: az(*) 
      real(double)  :: bx(*) 
      real(double)  :: by(*) 
      real(double)  :: bz(*) 
      real(double)  :: cx(*) 
      real(double)  :: cy(*) 
      real(double)  :: cz(*) 
      real(double)  :: dx1(*) 
      real(double)  :: dx2(*) 
      real(double)  :: dx3(*) 
      real(double)  :: tsix(*) 
      real(double)  :: tsiy(*) 
      real(double)  :: tsiz(*) 
      real(double)  :: etax(*) 
      real(double)  :: etay(*) 
      real(double)  :: etaz(*) 
      real(double)  :: nux(*) 
      real(double)  :: nuy(*) 
      real(double)  :: nuz(*) 
      real(double) , intent(out) :: divaa1(*) 
      real(double) , intent(out) :: divbb2(*) 
      real(double) , intent(out) :: divcc3(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: istart, jstart, ncells, i, j, k, ijk, iter, n 
      real(double) :: jacob 
!-----------------------------------------------
!
!
!
      istart = 3 
      jstart = 3 
      if (periodic) then 
         istart = 2 
         jstart = 2 
      endif 
!
!     construct an array of cell indices
!
      ncells = 0 
      do i = istart, nxp - 1 
         do j = jstart, nyp - 1 
            do k = 3, nzp - 1 
!
               ncells = ncells + 1 
!
               ijkcell(ncells) = (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid + 1 
!
               ijk = ijkcell(ncells) 
               divaa1(ijk) = 0.0 
               divbb2(ijk) = 0.0 
               divcc3(ijk) = 0.0 
!
            end do 
         end do 
      end do 
!
!
      do iter = 1, numit 
!
         call contras (iwid, jwid, kwid, ncells, ijkcell, x, y, z, tsix, tsiy, &
            tsiz, etax, etay, etaz, nux, nuy, nuz) 
!
!     compute weight function for toroidal mesh
!
         call torusw (nxp, nyp, nzp, 1, nxp, 1, nyp, 2, nzp, toroid, strait, &
            rmaj, x, y, z, w) 
!
         call adaptw (nxp, nyp, nzp, 1, nxp, 1, nyp, 2, nzp, toroid, strait, &
            rmaj, x, y, z, w, w_soln) 
!
!
!     calculate change in covariant components of mesh vector
!
         call gridsolv (x, y, z, w, tsix, tsiy, tsiz, etax, etay, etaz, nux, &
            nuy, nuz, iwid, jwid, kwid, ijkcell, ncells, dx1, dx2, dx3) 
!
!dir$ ivdep
         do n = 1, ncells 
!
            ijk = ijkcell(n) 
!
            x(ijk) = x(ijk) + 0.5*dx1(ijk) 
            y(ijk) = y(ijk) + 0.5*dx2(ijk) 
            z(ijk) = z(ijk) + 0.5*dx3(ijk) 
!
         end do 
!
!     impose periodic boundary conditions for torus
!
         call torusbg (nxp, nyp, nzp, istep, iota, cdlt, sdlt, strait, toroid, &
            rwall, rmaj, dz, x, y, z) 
!
      end do 
!
      return  
      end subroutine meshgen 


 
      subroutine torusw(nxp, nyp, nzp, i1, i2, j1, j2, k1, k2, toroid, strait, &
         rmaj, x, y, z, wgrid) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp 
      integer , intent(in) :: nyp 
      integer  :: nzp 
      integer , intent(in) :: i1 
      integer , intent(in) :: i2 
      integer , intent(in) :: j1 
      integer , intent(in) :: j2 
      integer , intent(in) :: k1 
      integer , intent(in) :: k2 
      real(double) , intent(in) :: toroid 
      real(double) , intent(in) :: strait 
      real(double) , intent(in) :: rmaj 
      real(double) , intent(in) :: x(nxp,nyp,*) 
      real(double) , intent(in) :: y(nxp,nyp,*) 
      real(double) , intent(in) :: z(nxp,nyp,*) 
      real(double) , intent(out) :: wgrid(nxp,nyp,*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, j, i 
      real(double) :: ws, wsr, wsmr 
!-----------------------------------------------
!
      do k = k1, k2 
         do j = j1, j2 
            do i = i1, i2 
!
               ws = toroid*sqrt(x(i,j,k)**2+y(i,j,k)**2) - rmaj + strait*x(i,j,&
                  k) 
               wsr = sqrt(ws**2 + z(i,j,k)**2) + 1.E-20 
               wsmr = toroid*sqrt(x(i,j,k)**2+y(i,j,k)**2) + strait 
!
!     to prevent collapse near the origin
!
!      wsrp=1./(1./wsr+0.25/wsr**2)
!
               wgrid(i,j,k) = wsr*wsmr 
!
            end do 
         end do 
      end do 
!
      return  
      end subroutine torusw 
